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INTRODUCTION 



I . 



A. BACKGROUND 

Mankind has used the ideas of science and mathematics to 
help understand, describe and control the world. The central 
notions of science have developed gradually. Truly revolu- 
tionary ideas do not often surface. The development of 
calculus by Newton and Leibnitz in the seventeenth century is 
quite possibly the last great revolutionary development in the 
sciences. Chaos and chaotic dynamics could very well be the 
next scientific revolution. The theories of chaotic dynamics 
are being used to describe and interpret many natural phenome- 
non that have heretofore been inexplicable. 

The modern study of chaotic dynamical systems was begun by 
Henri Poincare at the beginning of the twentieth century. 
During his study of celestial bodies, Poincare observed that 
their trajectories were forever oscillating, yet irregular and 
aperiodic. Poincare described this motion as chaotic. He 
further discovered that the motion of a chaotic system has a 
very sensitive dependence on the initial conditions of the 
system. Because of the nonlinearity of chaotic systems, 
further rigorous study was impossible until the advent of the 
computer. Since 1960, many mathematicians and engineers have 
developed techniques to analyze chaotic systems. Much debate 
has occurred about the nature of chaos. Currently, the only 
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universally accepted property of chaotic dynamical systems is 
that they have sensitive dependence to initial conditions. 
Agreement about certain aspects of chaos does not exist within 
either the mathematical community or the engineering community 
and certainly no agreement exists between these groups. 

In its David II report [Ref 1.1], the Committee on the 
Mathematical Sciences of the National Research Council, listed 
chaotic dynamics as one of the brightest areas for future ad- 
vancement of applied mathematics. The report states that 
chaotic behavior has been found in nonlinear systems as 
diverse as ecology and economics to engineering and meteorolo- 
gy. The theories of chaos have their basis in the mathemati- 
cal fields of topology, measure theory and combinatorics. 
Mathematicians are just now beginning to develop numerical 
methods for computing some of the properties which Poincare 
described some eighty years ago. 

The engineering community looks at chaos from a different 
viewpoint. Engineers have developed techniques to analyze 
systems that have discrete data points to determine their 
chaotic tendencies. Pseudo phase planes, Lyapunov exponents 
and multivariate probability scaling analysis are but a few of 
the techniques that lend insight into the behavior of chaotic 
dynamical systems. No single method is completely conclusive. 
Hence, when studying nonlinear systems, a variety of methods 
is used. 
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The motivation for this thesis was twofold. First, an 
attempt has been made to bridge the gap between how mathemati- 
cians develop chaotic analysis and how engineers apply the 
methods of chaotic analysis. This was accomplished by 
examining systems that mathematicians have proven to be 
chaotic and using engineering methods to analyze these 
systems. The fortran program CHAOS developed by 
CDR Martinus Sarigul-Kli jn [Ref. 1.2] while studying heli- 
copter vibrations was used to analyze the Lorenz equations and 
Duffing' s equation. Secondly, the data generated from the 
Lorenz equations. Duffing' s equation and the helicopter 
vibration data used by CDR Sarigul-Kli jn was analyzed using a 
new technique proposed by Osborne and Provenzale [Ref 1.3]. 
All computations were performed on a VAX based system in the 
Advanced Computational Laboratory of the Department of 
Aeronautics and Astronautics, Naval Postgraduate School. 

B. OVERVIEW OF THESIS 

The second chapter of the thesis covers background needed 
to understand chaos. The chapter begins with a short develop- 
ment of the mathematical basis for chaos. Probability theory, 
specifically the central limit theorem, is also developed. 
Then, the engineering techniques used to analyze chaos are 
presented. Pseudo phase planes, fractal dimension, Poincare' 
sections and Lyapunov exponents are presented here. Finally, 
the chapter ends with a literature review that discusses 
recent developments and relevant papers concerning 
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applications of chaos. 

The third chapter begins with a discussion of a technique 
for simulating a continuous system using a digital computer. 
Next, the work of CDR Sarigul-Kli jn is presented. Next, the 
methods used in the present analysis are presented. The 
multivariate scaling analysis is presented here. Further, the 
method used to analyze a continuous system using a discrete 
approximation is discussed. 

Chapter IV presents a summary of the results of the 
analysis. First, the Lorenz equations and Duffing' s equation 
are evaluated using current techniques. Next, the multivari- 
ate scaling analysis results are presented for the helicopter 
vibrations data. 

Chapter V is a discussion of conclusions from the analy- 
sis. Scope for further research is also presented. 
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UNDERSTANDING CHAOS 



II . 



A. INTRODUCTION 

This chapter gives a brief overview of the background 
needed to understand the concept of chaos. First, the 
mathematical basis for choatic dynamical systems is presented. 
Next, the engineering techniques used to analyze chaos are 
presented. Finally, the chapter ends with a review of recent 
literature that concerns applications of chaos that are 
relevant to the present study. 

B. MATHEMATICAL BASIS FOR CHAOS 

The study of chaotic systems begins with the study of dis- 
crete dynamical systems. Since data can only be collected at 
discrete points in time, it is quite natural to examine a 
dynamical system at those discrete points so that a model of 
the system can be constructed. The computer has allowed the 
analysis to be done with shorter and shorter time intervals. 
A simple example of a discrete dynamical system is an iterated 
function system constructed with a calculator. A value in 
this system depends on the previous value in the system. The 
system could be constructed by choosing a value and then 
repeatedly pressing the same function key on a calculator. If 
the key selected were the square root function and the 
starting number was 100, then we could construct a sequence of 
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numbers based on the function (square root) and the initial 
condition. Clearly this system is quite predictable but it is 
still dependent on the initial conditions. A system with this 
property is said to be deterministic. 

Dynamical systems can be classified by the nature of the 
model that best describes that system. A 'linear' system is 
modeled by a funtion that is a straight line that passes 
through origin. A function is called 'affine' if its graph is 
a straight line that does not pass through the origin. An 
affine system could be modeled by the relationship: 

= AX„ + B. 

A 'nonlinear system' takes the form: 

X„,i=f (XJ+g(n) 

where f (X^) = AX„ and g(n) is a function of the iteration 
number. The graph of a. nonlinear system is not a straight 
line and need not pass through the origin. For example, if 
f (X) = X and g(n) = n the equation for the 15th iteration 
would be 

Xi5=f (Xi4)=Xi4 +15. 

A simple iterated function system might revisit a point 
after a certain interval. If X^+i = f (X„) = X„ then X„ is 
called a 'fixed point' . If X„ = f (X^) = X„+„ then X„ is called 
a 'periodic point'. That is, if one starts at point X„ and 
performs 'm' iterations on the system the last iteration will 
yield point X„. A fixed point remains fixed after successive 
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iterations, while a periodic point is revisited only after a 
certain number of function iterations. 

Periodic points and fixed points can be classified as 
attracting or repelling. A point is called 'repelling' if a 
small perturbation in the initial condition causes the 
function to diverge exponentially. A point X is called at- 
tracting if a small perturbation in the initial condition 
causes the system to converge to point X after successive 
iterations. That is, if a system is iterated with two sets of 
initial conditions, one condition an attracting point and one 
condition e distant from the attracting point, then after 
successive iterations the functional values will be closer 
together than the original points. Figure 2.1a and 2.1b 
exhibit a repelling point and an attracting point, respective- 
ly. It can be shown that if | f ' (X^) | < 1 then X„ is an attract- 
ing point. While, if | f ' (X„) | > 1 then X„ is a repelling 
point . 

Before giving a somewhat generally accepted mathematical 
definition of chaos, the concepts of transitivity, denseness 
and other basic mathematical concepts must be explored. Many 
good references are available. Devaney [Ref. 2.1] presents a 
particularly good explanation of these mathematical concepts. 

First the definition of metric space is presented. A 
'metric space' (X,d) is a space X, together with a real valued 
function d, that maps points in X to points in X, and which 
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measures the distance between pairs of points in X. The 
function d must obey the following rules: 

(1) d(x,y) = d(y,x) for all x and y in X 

(2) 0 < d(x,y) < oo for all x, y in X, x y 

(3) d(x,x) = 0 for all x in X 

(4) d(x,y) < d(x,z) + d(z,y) for all x, y, and z in X 

A point is called a 'limit point' if there is a sequence 
of points that converge to that point, A connection between 
attractors and limit points becomes evident. The concept of 
convergence is common to both attractors and limit points. 
Next, the 'closure' of a subset of a metric space is said to 
be the union of that subset and all of its limit points. 
Next, 'denseness' is defined. Let (X,d) be a metric space. 
Let B be a subset of X. B is said to be dense in X if the 
closure of B equals X. 

A set is 'closed' if it contains all of its limit points. 
If a set is not closed it is said to be 'open' . The concept 
of 'transitivity' is next presented. Let U and V be open 
subsets of the metric space (X,d) . A dynamical system {x,f} 
in that metric space is transitive if there exists a finite 
integer such that 

U C\ *0 

In simple terms, this means that any arbitrarily small 
neighborhood in X will, under successive iterations, overlap 
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with any other arbitrarily small neighborhood in X. There- 
fore, the system cannot be divided into two disjoint open 
sets . 

Now the definition of 'chaos' as defined in [Ref. 2.1] is 
presented. A system is said to be chaotic if: 

(1) it has sensitive dependence on initial conditions 

(2) it is topologically transitive 

(3) periodic points are dense. 

An important phenomena in the study of chaos is the 
'bifurcation' . In general, dynamical systems arrive at a 
chaotic state through a bifurcation of one type or another. 
Bifurcations are found throughout the real world. A bifurca- 
tion is simply a dramatic change in a dynamical system. Water 
freezing, beam buckling under loads and the splitting of an 
atom during a nuclear explosion are all examples of bifurca- 
tions. Many systems are nonchaotic until a parameter or set 
of parameters attain certain critical values. As the system 
achieves a critical state a bifurcation occurs and the system 
becomes chaotic. In general, not all bifurcations are identi- 
cal . 

A 'saddle-node bifurcation' occurs when a function has a 
semi-stable fixed point. A semi-stable fixed point is 
attracting if perturbed in one direction and repelling if 
perturbed in another direction. In the simplest two dimen- 
sional case, the bifurcation occurs as the function approaches 
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the line y =x. When the graph of the function just touches 
the line the function splits into two. Figures 2,2a and 2.2b 
provide graphical representations of a saddle-node bifurcation 
of the logistic equation. 

A period doubling bifurcation occurs when an attracting 
fixed point becomes repelling. As the parameter achieves its 
critical value, the fixed point now becomes a cycle of period 
two. As the parameter passes another critical value, yet 
another split occurs. This process continues until a chaotic 
state is achieved. Saddle-node and period doubling are the 
two most common bifurcations. While others are found in 
nature, it is clear that a route to chaos is through 
bifurcations [Ref. 2.2]. 

A remarkable property about dynamical systems was proved 
by Sarkovskii [Ref 2.1]. Sarkovskii developed a procedure to 
determine the periods of a dynamical system. He showed that 
by ordering the integers in a certain manner, the periodic 
points of a function could be predicted. As previously 
discussed, a periodic point is a point that is revisited after 
a number of functions iterations. In this ordering if a 
system has a periodic point of a higher order period then it 
has periodic points of all the periods that follow. For 
example, if a system has a point of period three, then it has 
periodic points of all periods. To find points of other 
periods one must change the initial conditions. The highest 
order period is three and the lowest order period is one. 
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Sarkovskii's theorem states as follows: 



let the ordering of the integers representing periodic 
points of a system be: 

3>5>7>9 

2x3>2x5>2x7> 

2^x3>2^x5> 

2 ®> 2 ^> 2 ^> 22 > 2>1 

then, if a function, f, has a periodic point of period n, and 
n > K in the ordering, then f also has a periodic point of 
period K. Hence, period 3 implies all other periods. Period 
2 only implies period 1. 

The next concept to be discussed in this section, is the 
notion of a fractal. Fractals are subsets of metric spaces. 
Fractals have become important in helping to classify and 
quantify chaos. Fractal dimension is a useful tool in 
measuring chaos.' 

The mathematics behind the fractal is quite rigor- 
ous. [Ref. 2. 3] A basic example will present a clear picture of 
a fractal. The Cantor middle thirds set is the archetypical 
fractal . This set is constructed by starting with the 
interval [0,1] . The function removes the middle third of each 
remaining segment. After successive iterations the fractal 
remains. Figure 2.3 illustrates the Cantor set. What remains 
is the attractor for the system. Each segment is similar to 
the previous segment. This important property is called self- 
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similarity. The function for generating a fractal can be a 
rotation, a translation, a reflection or any combination of 
those three. Many functions or mappings may be needed to 
generate a particular fractal. In all cases, the function (or 
functions) acts on the previous point or segment. However, to 
generate a fractal each function must be contractive rather 
than expansive. This means that the function cannot enlarge 
the previous segment. But, the combination of multiple 
functions can be larger than the original segment. 

'Fractal dimension' provides a useful measure of a chaotic 
system. A chaotic system will have a noninteger fractal 
dimension. While not all systems with noninteger fractal 
dimensions are chaotic, this is still a useful measure [Ref. 
1.3] . 

The fractal dimension is a measure of how much space an 
object occupies in its metric space. An e-ball about a point 
is simply all points less than e distance from the point. To 
cover a set, one should take the smallest number of e-balls 
that intersect all points in the set. The fractal dimension 
is defined as follows: 

Let A be an element of the Hausdorf space where (X, d) is 
the metric. For each e > 0 let N(A, e) denote the smallest 
number of closed balls of radius e > 0 needed to cover A. 

„ _ Lim LnjN(A,e) ) 
e-*0 Ln{l/e) 
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Then, element A is said to have fractal dimension D. The 
fractal dimension gives a measure of the relative size of the 
attractor of a chaotic system. 

The examples discussed above were all one dimensional in 
nature. However, the concepts of chaotic dynamics can easily 
be applied to systems of higher dimensions. 

The last concepts discussed in this section are the 
relevant areas of probability theory. The probability density 
function is the relative measure of the frequency of occur- 
rences. For a discrete random variable the summation of all 
values of the probability density function is unity. That is, 
if f (Xn) is the probability of occurrence of point x„ then 

N 

72=1 



This concept is intuitively obvious. The analog to summation 
for a continuous system is the integral; for a system with 
more than one random variable the concept does not change. 
For a multivariable system the summation over all random 
variables will equal unity. The summation over a single 
random variable of a multivariable system will always be less 
than or equal to unity. 

Data can be distributed in a variety of patterns. If each 
value in the range has an equal probability of occurring, the 
data is said to have a uniform distribution. If the data is 
symmetrically dispersed about a central point with more data 
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near the center of the range than at the ends of the range, 
the data has a Gaussian or normal distribution. Figure 2.5 
illustrates some common probability distributions. A key 
theorem allows data with any distribution to be transformed to 
normal distribution. The Central Limit Theorem follows: 

Theorem: Let X^, . . . X„ be independent and identically 
distributed random variables all with expectations = |i and 
variance equal to a 



then converges to a standard normal random variable. 
This theorem is applied by Osborne and Provenzale [Ref 1.3] 
when comparing chaotic systems with randomly generated 
systems . 

C . ENGINEERING TECHNIQUES 

Most naturally occurring phenomenon have a continuous 
nonlinear pattern. To study these events, engineers must take 
discrete readings over a time interval or across a distance. 
The very act of measurement applies a strain on the system and 
slightly changes the system. The process of discretizing the 
continuous data produces an imperfect model of the system. 
Again error is introduced into the data. This section will 
discuss the engineering techniques used in the study of 
chaotic dynamics. 
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1. Time Domain Analysis 

The three steps to convert an analog signal to a 
digital signal are sampling, quantizing and encoding [Refs. 
1.2 and 2.4]. Sampling a continuous signal produces a series 
of discrete values. A generally accepted practice is to 
sample the data at twice the rate of the uppermost frequency 
of the signal. The sampled data must be quantized. The 
spectrum is divided into different levels known as quantum 
levels. The signal is then compared with each quantum level 
and the reading at a specific time is recorded at the nearest 
quantum level . The maximum error introduced here is then one- 
half the quantum level. The next step is translating the 
quantum levels to the maximum number of levels available in 
the data recording equipment. This is know as machine 
precision. Again an error is introduced into the data. 

2. Frequency Domain Analysis 

An alternative method for analyzing a continuous data 
signal is a frequency domain analysis. This method uses a 
Fourier analysis to measure a function using harmonic compo- 
nents that have varying amplitudes, frequencies and phases. 
A signal in the frequency domain can be represented by two 
plots. A power spectral density graph is a plot that is often 
used to represent signals. The square of the amplitude is 
plotted against the frequency. 
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3. Phase Planes and Pseudo Phase Planes 



The concept of the phase plane was introduced by Poin- 
care. A phase plane represents a dynamical system by using 
two independent properties that best describe the system 
dynamics. At each point in time those properties are measured 
and then plotted on a two dimensional graph. At any point in 
time the system can be described by a single point on the 
graph. The resulting plot is used to analyze the system. If 
the system approaches a limit point or limit cycle over time 
it is said to have an attractor. The phase plane requires two 
separate measurements of a system. As was discussed earlier, 
each measurement introduces an error into the data. Also, 
collecting data for two variables will be more costly and 
might be physically impossible. Hence, a phase plane analysis 
may be impractical. 

Takens suggested a simpler technique for constructing a 
phase plane diagram [Ref. 2.5]. This technique is called a 
pseudo phase plane. A dynamical system is measured with only 
one variable. This procedure generates a single time series. 
A second time series is constructed by adding an embedding 
time to the original data. If the embedding time is two time 
periods, the third point in the original series becomes the 
first point in the new series. The data in the second series 
is then plotted versus the first series to construct a pseudo 
phase plane. The pseudo phase plane captures the system 
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dynamics without filtering the signal by integration or 
differentiation of the measured signals [Ref. 2.6]. Any 
dynamic property of the system can be used to represent the 
system. The choice of embedding values plays a critical role 
in the accuracy and interpretation of the analysis. A small 
embedding time produces a centralized data set while a large 
embedding time introduces error through edge effects. 
Grassberger [Ref. 2. 7] suggests that the embedding time should 
be about one quarter of the most predominant frequency of the 
data. Dvorak and Klaschka [Ref. 2.8] have since proven that 
an embedding time between five and eleven time intervals 
produces a pseudo phase plane with the least amount of bias. 
In recent years, many methods have been developed for 
reducing the measurement error or noise in the analysis. 
Hammel [Ref. 2.9] has been able to reduce the effects of noise 
tenfold. This produces a data set that is accurate to machine 
precision on even the most advanced computers. For an 
experimental procedure or a numerical simulation, the pseudo 
phase plane is more efficient, less costly and produces 
results that are just as accurate as the phase plane and 
captures the original system dynamics . 

Pseudo phase planes can be constructed in higher dimen- 
sions as well. The procedure to construct higher dimensional 
pseudo phase planes is similar to that used to construct the 
two dimensional phase planes. The third and subsequent time 
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series are constructed using the original time series and new 
embedding factors. The embedding factor should be different 
for each new time series. More than three dimensions cannot 
be depicted graphically but, multidimensional pseudo phase 
planes can provide insight into the dynamics of a system. 

Poincare sections can be constructed from these phase 
planes [Ref. 2.10]. A Poincare section is a two dimensional 
cut of a three dimensional phase plane. The section shows 
where trajectories of the system pierce the phase plane. This 
section allows one to graphically determine if a bounded 
attractor or limit cycle is present. A well defined attractor 
is evidence of a chaotic system. The Poincare section can be 
taken at any angle so that the number of possible sections is 
infinite. Higher dimensional Poincare sections may be 
obtained similarly. 

4 . Van Der Pol Planes 

In higher dimensional space, Poincare sections can be 
obtained as intersections of the trajectories with a specified 
hyper-plane. The internal structure of the attractor may then 
be examined. The Van Der Pol plot is constructed by untwisting 
the trajectory at a certain rate to create a series of 
Poincare sections. The sections are then merged together to 
yield a Van Der Pol plot [Ref. 2.11] . These are particularly 
useful when a single frequency dominates the response of a 
system. More information on Poincare sections can be found 
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in [Ref. 2.6] . 

5. Lyapunov Exponents 

The methods discussed previously have all provided a 
graphical representation of the chaotic nature of dynamical 
systems. The graphs offer pictorial evidence of the chaotic 
nature of the system being analyzed. These methods do not 
provide quantitative measures of chaos. The Lyapunov exponent 
offers quantitative evidence of chaos. Wolf [Ref. 2.12] 
offers an algorithm for calculating these exponents based on 
a time series. This method measures the attraction or repul- 
sion over time of two adjacent trajectories with different 
initial conditions. A positive exponent means that the 
distance between the systems has grown greater with time. A 
negative exponent indicates that the trajectories have grown 
closer. A positive exponent is an indication of unpredict- 
ability and chaotic behavior in a dynamical system. A 
positive exponent alone is not enough to determine if a system 
is chaotic. Recent studies by Kapitaniak and Naschie [Ref. 
2.13] show that a purely random system can also have positive 
Lyapunov exponents . 

The Lyapunov exponents also gives a measure of the rate of 
information loss in a system. As such, the exponent can be a 
useful gage as to how far into the future one can predict the 
dynamics of a system. 
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6. Fractal Dimension 



Another measure used to quantify chaos in a dynamical 
system is the fractal correlation dimension. As discussed 
earlier, fractal dimension gives a lower bound for the number 
of degrees of freedom needed to model a physical system. A 
finite non-integer fractal dimension suggests yet further 
evidence of the presence of a strange attractor and chaos. 
However, Osborne and Provenzale have shown that a randomly 
constructed data can also have a non-integer fractal dimension 
[Ref. 1.3]. Ruelle [Ref. 2.14] has shown that all fractal 
dimension estimates that approach 21ogN are suspect. In this 
case, N is the number of data points the time series. 
Therefore, for most time series of 10,000 data points fractal 
dimension approaching six yields inconclusive evidence of 
chaos in the dynamical system. 

As with much of this field, no universally accepted 
definition of fractal dimension exists as yet. Farmer et al 
[Ref. 2.15] describe seven natural measures of dimension that 
could be called fractal dimension. These measures can be 
divided into two broad classes. One class requires a proba- 
bility measure to be computed. The other class requires only 
a geometric measurement . 

7 . Evolving Measures of Chaos 

Osborne and Provenzale [Ref. 1.3] have suggested that 
a multivariate scaling analysis will lend additional insight 
into the dynamics of a system. A multivariate scaling 
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analysis involves computing probability density functions (PDF) 
of the time series data. The data is normalized before the 
PDF is constructed. A non-chaotic system will generate a 
Gaussian PDF. A chaotic system will yield a non-Gaussian PDF. 
These ideas will be further elaborated in Chapter III. 

Lindsay [Ref 2.16] suggests a method for forecasting data 
points of a chaotic system. This method calls for calculating 
the Jacobian and strange attractor to use global rather than 
local properties of the attractor. A predictor function is 
constructed using the Jacobian that allows one to forecast as 
many as one-half the original number of data points. 

The last two methods utilize the global properties of 
chaotic systems. Locally, random systems often mirror chaotic 
systems. A global approach helps to minimize the confusion 
caused when comparing local properties of random systems with 
local properties of chaotic systems. 

D. SUMMARY 

Techniques to analyze chaotic dynamical systems have 
emerged rapidly during the last decade. Mathematicians, 
physicists and engineers differ in their approach concerning 
analyzing the nature of chaotic systems. Methods to investi- 
gate chaos in systems modeled by differential equations are 
well documented. Such systems are continuous and need not be 
modeled with discrete data points. Measuring real physical 
systems yields discrete data points rather than continous 
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equations. The measurement process introduces several sources 
of error into the system. Many methods can help describe the 
chaotic nature of such a system. But, as a result of the 
error introduced, no method gives a fool- proof description of 
chaos. Many techniques must be employed concurrently to 
ascertain the evidence of chaos and its associated properties 
that allow accurate definition of the dynamical system. 

E . LITERATURE REVIEW 

Researchers in the field of chaotic dynamics are breaking 
new ground daily and technical papers concerning chaos 
proliferate in leading scientific journals. In order to limit 
the scope of this literature review, the papers and reports 
discussed in this section are those that were published 
between December of 1989 and April of 1991. This list is not 
exhaustive, but covers the papers that are directly related to 
the subject matter of this thesis. 

Bleher, Grebogi and Ott [Ref 2.17] have observed a new 
transition from a system in a steady-state to a system in a 
chaotic state. This transition is called an abrupt bifurca- 
tion. The abrupt bifurcation occurs in the context of chaotic 
scattering. The behavior is caused by the presence of a 
chaotic invariant set of unstable bounded orbits. An example 
showing the energy involved in atom-diatom collisions is 
discussed. Because of its dependence on a set of unstable 
bounded orbits, the abrupt bifurcation might help explain the 
nature of helicopter vibrations. 



22 



Hughes and Proctor [Ref. 2. 18] studied the effects of noise 
on a chaotic system. They show that noise, even at very low 
levels, can create great problems in determining the chaotic 
measures of a system. When noise was added to a known 
deterministic chaotic system the solution to the system is 
shown to be indeterminate and incomputable. 

Corless, Frank and Monroe [Ref. 2.19] determine the 
chaotic properties of the Gauss map. The Gauss map is 
developed using continued fractions. The authors determine 
the fractal dimension, Lyapunov exponent and the correlation 
dimension of the continous Gauss map. All of these properties 
indicate that the Gauss map is chaotic. Next, a numerical 
computer simulation is performed to generate a Gauss map. The 
simulated map has fractal dimension, Lyapunov exponent and 
correlation dimension that indicate the simulated system is 
nonchaotic. The simulated system is created to machine 
precision. The observable behavior of the systems is not 
significantly different. This suggests that a system with a 
long period might be chaotic. 

A numerical experiment using a self-affine time series is 
presented by Higuchi [Ref 2.20] . The randomness of the phase 
distribution strongly affects the irregularity of the system. 
This conclusion is based on a comparison of the phase distri- 
bution with the fractal dimension. The randomness of the 
phase distribution is shown to be affected by the behavior of 
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the system in the time domain. 

Liebert and Schuster [Ref. 2.21] report the effects of 
varying the embedding time delay when analyzing chaotic time 
series. For an infinite time series, any arbitrary value of 
an embedding value is acceptable when constructing a pseudo 
phase plane. However, for finite time series, arbitrary 
choices of embedding values can give erroneous results when 
calculating entropy and density. A method for calculating a 
proper embedding value which requires at least 9000 data 
points is presented. 

Fang [Ref. 2.22] suggests that the interplay between a 
large noise intensity and the nonlinearity of the system may 
induce a nonequilibrium transition of the dynamical behavior 
of the system. The degree of distortion is based on the ratio 
of the noise intensity to a measure of the system nonlinear- 
ity. The chaotic measures of the system may be affected by 
that ratio. Specifically, the Lyapunov exponent changes sign 
as the noise level is adjusted. The system changes from a 
chaotic state to a periodic state with the introduction of 
strong white noise. 

Kapitaniak [Ref. 2.23] presents a new method for estimat- 
ing the parameter values that cause a dynamical system to be 
chaotic. This method applies to systems that exhibit a period 
doubling route to chaos. The method is applied to Buffing's 
oscillator. The method is based on an approximate analysis 
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and the Feigenbaum universal properties of period doubling. 

Dvorak and Klaschka [Ref 2.8] present a modification to 
the Grassberger and Procaccia algorithm for finding correla- 
tion dimension. The original algorithm produces a bias for 
embedding dimensions greater than five. The correlation 
dimension of the system is changed because of the noise 
present in the algorithm. The new algorithm allows embedding 
dimensions up to eleven by reducing the white noise level. 
This algorithm removes the bias caused by the edge effects. 
All other causes of bias from the old algorithm are present in 
the new algorithm. 

Chen and Ott [Ref 2.24] present an algorithm for computing 
the cross-sections of strange attractors. These cross- 
sections are useful for elucidating the geometry of the 
attractor. The cross-sections also allow one to estimate the 
fractal dimension of higher dimensional attractors. 

Hammel [2.9] presents a method for reducing the noise of 
a time series of discrete data points. This method reduces 
the noise level tenfold. The reduction of noise allows one to 
distinguish the effects of noise from the effects of chaos. 
Both additive noise and dynamic noise levels are reduced. 
This method produces a time-series that is accurate to machine 
precision and thus yields more accurate calculations of 
chaotic measures of dynamical systems. 

Stone [Ref 2.25] compares the effects of noise versus 
periodic forcing on the power spectra of a Duffing oscillator. 
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It is shown that the form of the power spectral density of the 
noise driven system and the deterministic chaotic system are 
indistinguishable. This suggests that a high frequency power 
spectrum cannot distinguish between a deterministic and 
stochastic system. 

Lindsay [Ref. 2.16] presents a simple, computationally 
efficient method for forecasting chaotic time series using the 
Jacobian and certain weighting functions. The method requires 
solving a small number of linear equations. If the original 
time series has 10,000 data points, this method can calculate 
the next 5000 data points with exceptional accuracy. The 
error for the 5000th point is shown to be less than 2.7x10"^%. 
The same method can be used for noise reduction. 
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Figure 2.1b Attracting Point 
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Figure 2.2 Period doubling route to chaos: 
logistic equation F (a) =r (a) (l~a) 
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Figure 2 . 3 
A Typical Fractal 
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Ill, METHODOLOGY 



A. INTRODUCTION 

This chapter presents the methods used to generate the 
simulated data and the methods used to study the nonlinear 
dynamical system data simulated and one set of helicopter 
flight vibration data. The first section presents the 
technique used to generate a discrete representation of a 
continuous system for the Lorenz equations and Duffing' s 
equation. A discussion of randomly generated time series is 
also presented. Next, Sarigul-Kli jn' s work is discussed. 
This section focuses on the calculation of the fractal 
dimension and the Lyapunov exponent of a dynamical system 
represented by a time series. A development of the multi- 
variate scaling analysis is then presented. The chapter 
closes with a review of the computer simulations used in the 
current analysis . 

B. SIMULATING A CONTINUOUS SYSTEM 

In 1963, Lorenz with the help of a computer, studied the 
atmosphere of a convecting fluid modeled by the now famous 
Lorenz equations . [Ref . 2.14] These equations were used 
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dx^ 

~dt 



=a 



(X2-X1) 



dx, 

— =pXi-X2-XiX3 



dx. 



3 _ 



dt 



= -jSx^ 



where 

a is the Prantl number 
p is the Rayleigh number 
fi is the aspect ratio. 

Lorenz studied these equations with the following combination 
of parameters 
G = 10 
p = 28 
li = 8/3 

used to model the system. Lorenz argued that the atmosphere 
modeled by this system of equations is chaotic and therefore, 
with current methods longterm behavior could not be predicted. 
The earth's atmosphere is similar to the system modeled by the 
Lorenz' equations and thus, longterm weather patterns are also 
unpredictable. Lorenz used the term chaotic to mean that the 
system had sensitive dependence on initial conditions. In 
1980, Ueda [Ref. 3.1] studied the Duffing equation to model 
the motion of a system with a forcing function undergoing 
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large elastic deflection. Ueda found that the time series 
generated from the behavior of this equation 

+0 . 05 +x^=7 . 5cos ( t) 

dt2 at 

also exhibited chaotic tendencies. That is, a long term 
prediction about the behavior of the system is impossible. 
Both buffing's equation and the Lorenz equations are models 
of continuous dynamical systems. These equations also model 
the dynamics of several systems. 

To effectively analyze these systems with a digital 
computer a time series of discrete data points is needed. In 
the present study, a fourth order Runge-Kutta numerical 
integration was used to generate the time series. Many 
excellent references are available for this and other numeri- 
cal integration schemes [Refs. 3.2 and 3.3] . The fourth order 
Runge-Kutta technique was used because it provided the level 
of precision required for the analysis. An adaptive step size 
control procedure was not used because the adaptive strategy 
is not compatible with the use of an embedding technique for 
generating a pseudo phase plane. To accurately generate a 
pseudo phase plane the original time series should have 
uniformly spaced data. Adaptive step size does not provide 
evenly spaced data and therefore was impractical for the 
current analysis. Acceptable precision was achieved by using 
a step size of .01 seconds. The equation parameters, step 
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size, embedding dimension and initial conditions were all 
varied to generate different time series for analysis. The 
time series were then analyzed. 

Normally distributed colored random noise was numerically 
simulated using a random number generator. The generator 
yielded uniformly distributed random numbers. The Central 
Limit Theorem was applied to those numbers to generate a 
series of normally distributed numbers. Twelve uniform random 
numbers were used to produce one Gaussian random number [Ref. 
3.4] . The numbers were then added to a sine function to 
produce a time series of simulated colored noise. The 
equation used to generate the series is 

X= . 9 +. 25sin {2n) + .IRAND 

C. AN OVERVIEW OF SARIGUL-KLIJN' S WORK 

Sarigul-Kli jn performed an analysis of helicopter flight 
vibration data [Ref. 1.2]. This section is a brief descrip- 
tion of part of his work. Although helicopter vibrations are 
mostly periodic, evidence of chaos was claimed. Specifically, 
the Lyapunov exponent for the vertical acceleration measured 
under the pilot's seat was found to be 0.3 to 1.7 bits per 
second. Also, the fractal dimension for this data was found 
to converge to 6.6 as the embedding dimension was increased. 
The fractal dimension for the lateral and longitudinal 
accelerations converged to 6.3 and 6.4, respectively. A 
random signal was found to have a linearly increasing fractal 
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dimension. This random signal had uniform distribution, 
instead of the more widely used normal distribution. 

The Lyapunov exponents were calculated based on an algori- 
thm developed by Wolf et al. [Ref. 2.12]. The exponent is 
defined as follows 



A=log2 




where 

d(0) is the initial distance between trajectories 
d(t) is the distance at time ' t' 

The algorithm measures the distance between points in the 
original time series and the points in the embedded time 
series. A positive exponent for any measurable observable 
indicates diverging separation between trajectories and hence 
evidence of chaos. 

Sarigul-Kli jn adopted the method attributed to Grassberger 
and Procaccia to calculate fractal correlation. This is one 
of many measures of fractal dimension [Ref. 2.15]. This 
measure of fractal dimension is defined as follows: 

C(r) 
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where 



C(r) is the probability of the attractor within a ball 



of radius r (r-ball) 



r is the radius of the ball 



d is the fractal dimension. 



Therefore 




In C(r) 
In r 



The fractal correlation dimension is computed by dividing the 
number of points inside the r-ball by the total number of 
points in the attractor. The process is repeated for several 
values of ' r' .The slope of the In C(r) versus In r curve gives 
the value of the fractal correlation dimension. 

The fractal dimension is obtained by applying this 
procedure to all dimensions of the pseudo phase space. As the 
phase space dimensions increase, the correlation dimension 
will approach an asymptotic value. For chaotic systems, this 
value should be a finite noninteger value. 

D. MULTIVARIATE SCALING ANALYSIS 

This section presents the multivariate scaling analysis 
proposed by Osborne and Provenzale [Ref. 1.3]. Osborne and 
Provenzale found that a stochastic process, simple colored 
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random noise, can have a finite fractal correlation dimension. 
Previously, a finite fractal correlation dimension was thought 
to be positive indication of the presence of chaos. A new 
method, the multivariate scaling analysis, was developed to 
aid in the search for chaos. 

The analysis begins with a time series of data. A second 
time series is generated with an embedding procedure or by 
taking a second measurement of the same variable. In the 
present analysis, a time embedding procedure was used to 
generate the second time series. The two time series are used 
to construct a scaling space. The coordinates for the scaling 
space are generated as follows 

. for i*j (12) 



where 

Sn are the scaling space coordinates 
Xn(ti) is the original time series 
Xn(tj) is the embedded series 

M is the number of points in the original series. 
Since the points where i^^j are not used, the time series of 
scales will have data points. The embedding dimension 

determines how many sets of scaling coordinates are generated. 
The dimension of the scaling space is equal to the embedding 
dimension. If the embedding dimension is two, then two sets 
of scaling coordinates are generated. A second series of 
scaling coordinates is constructed using the same procedure 
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but a different embedding value than the first series. A 
joint probability density function is then constructed using 
the procedure discussed in Chapter II. This process can be 
applied to create a scaling space of two or higher dimensions. 
However, a graphic analysis only allows two dimensions to be 
viewed. The probability density function is analyzed to 
determine the chaotic nature of the system. A nonchaotic 
system will yield a joint probability density function that is 
Gaussian or nearly Gaussian in nature. A chaotic system will 
yield a joint probability density function that is clearly 
nonGaussian in nature. 

In the present analysis, multivariate scaling spaces and 
joint probability density functions were constucted for the 
Lorenz equations, buffing's equation, helicopter vibration 
data and colored random noise. While the multivariate scaling 
analysis does not provide conclusive proof of the presence of 
chaos, it does provide evidence to support a thorough investi- 
gation of a dynamical system. 

E . COMPUTER SIMULATIONS 

The analysis of data was conducted using the methods de- 
scribed in previous sections. Time series were numerically 
generated for the Lorenz equations, buffing's equation and 
colored random noise. Many different parameter values were 
studied for the Lorenz and buffing's systems. bifferent 
values of step size were tried for the numerical integrations. 
For all cases, effect of changes in embedding value was 
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studied. Lastly, for all cases, various sets of initial 
conditions were evaluated. 

The evaluations consisted of observing time history, 
computing fractal correlation dimension, Lyapunov exponents, 
Fourier spectrum and conducting a multivariate scaling space 
analysis for each set of parameters, embedding values, initial 
conditions and step size. Each variable could conceivably 
affect the chaotic nature of the system. Significant results 
are presented in Chapter IV. 
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IV. RESULTS 



A. INTRODUCTION 

This chapter presents the results of the analysis of the 
Lorenz equations, Duffing' s equations, random colored noise 
and helicopter flight vibration data. The results are first 
presented for each dynamical system. In order to study the 
dynamical system described by a system of ordinary 
differential equations, a program was developed based on a 4th 
order Runge-Kutta integration scheme. The program generated 
the time dependent information necessary to study the 
dynamical system. Another program, based on the multivariate 
scaling analysis, uses the time histories to compute joint 
probability density functions as discussed in earlier 
chapters. Listings of these programs are included in 
Appendices A and B, respectively. 



B. RANDOM COLORED NOISE 



Random colored noise is generated using the following 
equation : 



X= . 9 + .25 sin {2n) +.1RAND 



The data generated is normalized to yield a Gaussian 
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distribution. A comparison of fractal correlation dimension 
for random colored noise and purely random data is included as 
Figure 4.1. The number of points generated is 10,000 with an 
embedding time of 10 samples. The fractal correlation 
dimension of the random data increases linearly to infinity 
with the embedding dimension. The fractal correlation 
dimension of the random colored noise increases at a slower 
rate than the embedding dimension and appears to reach an 
asymptotic limit. The fractal dimension of random colored 
noise is a noninteger finite number for all values of 
embedding dimension that were tested. This data supports the 
conclusion of Osborne and Provenzale [Ref. 1.3] . A multivari- 
ate scaling analysis of random colored noise yields a joint 
probability density function that is seen to be Gaussian in 
nature. A graph of this function is included as Figure 4.2. 

C. LORENZ EQUATIONS 

The time series for the Lorenz equations is numerically 
generated using a fourth order Runge-Kutta integration 
numerical scheme. Step sizes of .001 seconds, .01 seconds and 
0.1 seconds were used. This system of equations could not be 
solved, even using double precision arithmetic, with a step- 
size of greater than .115 seconds on either the Micro-VAX II 
computer at the Naval Postgraduate School or the CRAY-YMP 
computer at NASA-Ames Center. Step sizes larger than .115 
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caused an arithmetic overflow error. Hence, all computations 
were performed with step-sizes less than .115 seconds. Figure 
4.3 shows a comparison of the X-Z phase plane with a pseudo 
phase plane generated using the X time series only and an 
embedding value of one. 10, 000 data points were used to 
generate the pseudo phase plane. The pseudo phase plane has 
the same general structure as the phase plane with a rotation 
of one branch of the diagram. Larger embedding dimensions 
caused the psuedo phase plane to rotate and generate more 
branches. The fractal correlation dimension increases with 
step size. Figure 4.4 shows a comparison of fractal dimension 
to step size. A step size of .1 seconds yields a fractal 
correlation dimension of 2.04 which agrees with the work of 
Grassberger and Procaccia [Ref. 4.1]. Figure 4.5 shows plots 
of the slope of the log(r) versus log C(r) curve with 
differing step sizes. 

An interesting result was observed as the Lorenz equations 
were studied. Smaller step sizes resulted in a loss of low 
frequency cycles in the Fourier spectrum. Figure 4.6 shows 
the Fourier spectrum plotted to the Nyquist frequency for the 
differing step sizes. Figure 4.7 shows plots to 20% of the 
Nyquist frequencies. Figure 4.8 shows the low frequency 
Fourier spectrum for different step sizes. The smaller step 
size appears to act as a filter for the low frequency 
components of the Fourier spectrum. The multivariate scaling 
analysis, using a 10,000 data points and a time delay of one 
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sample, exhibited a non Gaussian behavior suggesting a chaotic 
attractor [Ref 1.3] . The probability distribution function is 
included as Figure 4.9. 

D. DUFFING' S EQUATIONS 

Next Duffing' s equation, which is known to exhibit chaotic 
behavior is studied. Numerical solution was obtained by 
adapting the program developed for the Lorenz equations. 
Figure 4.10 shows a comparison of the X-Y phase plane versus 
the pseudo phase plane generated with an embedding dimension 
of 1 applied to the X time series. Again, fractal correlation 
dimension changed with step size. However, in this case the 
fractal dimension decreased with larger step sizes. In all 
cases, the fractal correlation dimensions were finite non- 
integer values. The range of fractal dimensions agrees with 
the work of Moon and Li’ [Ref. 4.2] . Figure 4.11 shows a 
comparison of fractal dimension versus change in step size 
using 10, 000 data points with an embedding dimension of 10 and 
a time delay of 1 sample. Figure 4.12 shows the Fourier 
spectrum plotted to the Nyquist frequency for the various step 
sizes. Figure 4.13 shows the Fourier spectrum plotted to 20% 
of the Nyquist frequency. Figure 4.14 shows a comparison of 
step size versus the low frequency Fourier spectrum. Again, 
the smaller step size acts as a low pass filter. But in this 
analysis the larger step size appears to cause the data to 
have quasiperiodic behavior. Figure 4.15 shows a graph of the 
joint probability density function generated from the 
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mut livariate scaling analysis. The graph is seen to be non- 
Gaussian indicating chaotic behavior. 

E. HELICOPTER VIBRATION DATA 

This section presents an application of the multivariate 
scaling analysis based on the pseudo phase plane method for 
the helicopter vibration data [Ref. 1.2]. The vibration data 
selected for the present study are the vertical acceleration 
recorded under the right pilot seat of an OH-6A helicopter. 
Both Higher Harmonic Control 'off' and 'on' are studied. The 
joint probability density function with a time delay of 10 
data points is included as figure 4.16. The graph differs 
little from a standard Gaussian distribution. The variation 
is more like the random colored noise than Lorenz or Buffing's 
equations. Figure 4.17 shows the joint probability density 
functions for all of the dynamical systems studied. This 
suggests that part of the helicopter vibrations are periodic 
while part of the vibration is chaotic. This conclusion 
supports the work of Sarigul-Kli jn [Ref 1.2]. 

F . CONCLUSIONS 

Two trends are evident in the analysis of the data. It 
appears that integration step size filters the data generated 
by numerical integration using a computer. Small step sizes 
might cause a loss of the low frequency components while large 
step sizes could cause the behavior of a chaotic dynamical 
system to appear to be quasiperiodic . The multivariate 
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scaling analysis appears to give an indication of the chaotic 
nature of a dynamical system distinguished from random data. 
Still, in the search for chaos, this technique probably does 
not provide a definitive answer as to the nature of chaotic 
behavior . 
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Figure 4.1 Fractal dimension vs. embedding 
dimension for random data and random colored 
noise . 
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Figure 4.2 Joint probability density 
function for 2-D scaling space generated 
from random colored noise using a multi- 
variate scaling analysis. 
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Figure 4.3a X-Z phase plane generated 
from Lorenz equatlofis. 
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Figure 4.3b X-x p.seudo phase plane 
generated from equations. 
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Figure 4,4 Fractal dimension vs. 
step size 
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Figure 4.5 Slope of log(r) vs. 
log C(r) for Lorenz equations. 
Maximum point is fractal correlation 
dimension . 
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Figure 4.6 Fourier spectrum for the 
Lorenz equations plotted to the Nyquist 
frequencies . 
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Figure 4.7 Fourier spectrum for the 
Lorenz equations plotted to 20% of 
the Nyquist frequencies. 
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Figure 4.8 Low frequency Fourier 
spectrum generated from Lorenz equations. 
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Figure 4.9 Scaling space joint 
probability density function for 
the Lorenz equations. 
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Figure 4.10n X-Y phase plane for 
Puffing's equation. 




Figure 4.10b X-X pseudo phase plane 
for Puffing's equation. 
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Figure 4.11 Slope of log(r) vs. log C(r) 
for Puffing's equation. The maximum 
point on the curve is the fractal 
correlation dimension. 
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Figure 4.12 Fourier spectrum plotted 
to the Nyquist frequency for 
Duffing' s equation. 
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Figure 4.13 Fourier spectrum for 
Duffing' s equation plotted to 
20% of the Nyquist frequency 
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Figure 4.14 Low frequency Fourier 
spectrum for Duffing' s equation. 
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Figure 4.15 Joint probability 
density function for scaling 
space for Duffing' s equations. 
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Figure 4.16 Joint probability 
density function for scaling 
space for helicopter vibration 
vertical acceleration data. 
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Figure 4.17 Scaling space joint 
probability density functions 
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V. CONCLUSIONS AND SCOPE FOR FUTURE RESEARCH 



Engineering methods of chaos theory were applied to a 
discrete representation of chaotic systems of equations. The 
Lorenz equations and Duffing' s equation were examined using a 
time series approach. The time series were generated numeri- 
cally using a fourth order Runge-Kutta integration scheme. 
Several simulations were executed varying the equation 
parameters, the initial conditions, the embedding dimension 
and the integration step size. Adaptive step size control was 
not used because it is not compatible with the use of a time 
embedding procedure to generate a pseudo phase plane. 
However, time stepsize affected the Fourier spectrums of the 
systems analyzed. Small step size acted as a low pass filter 
in the chaotic dynamical systems. For certain parameters, the 
systems produced Lyapunov exponents and fractal dimension that 
indicated the presence of chaos. The size of the embedding 
dimension was found to affect the chaotic nature of the 
systems. The best embedding dimensions were found to be 
between one and ten. Large embedding dimensions caused the 
pseudo phase plane to have a different nature than the phase 
plane . 

A new method, the multivariate scaling analysis, was 
applied using embedding coordinates to the Lorenz equations 
and buffing's equation as well as to the helicopter vibration 
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data studied by Sarigul-Kli jn . A uniformly distributed random 
time series and a normally distributed random time series were 
also examined using this new technique. When examined using 
the multivariate scaling analysis, chaotic systems yielded 
probability density functions that were nearly Gaussian. The 
normally distributed random time series produced a density 
function that was Gaussian in nature. Lastly, the helicopter 
vibration data yielded a density function that was non- 
Gaussian, indicating a small amount of chaos in helicopter 
vibrations . 

Many avenues are open to further research in this area. 
Techniques have been developed that allow accurate prediction 
of future points in a chaotic time series. This technique can 
be used to predict the next 5000 points in a time series for 
helicopter vibrations. This technique could be applied to the 
currently available data. A higher harmonic control device 
utilizing this technique could be used to actively dampen the 
helicopter vibrations. The multivariate scaling analysis 
could be used to examine other experimental data that appears 
to have a chaotic nature. This technique offers further 
evidence of the chaotic nature of the system. 

The Naval Air Test Center will soon instrument an H-60 
helicopter with a higher harmonic control system. The methods 
described herein could be applied to the data obtained from 
that study. This would offer further evidence of the chaotic 
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nature of helicopter vibrations and help to determine if 
vibration data can be accurately predicted and hence signifi- 
cantly reduced. 
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APPENDIX A 

PROGRAM TO PERFORM NUMERICAL INTEGRATION 



A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 



This program solves systems of ordinary differential equations 
using a fourth order Runge-Kutta Scheme. 

The program was used to solve both Duffing's equation and the 
Lorenz equations. 



Two subroutines are employed. 

DERIVS calculates the values for the ODE ' s 
RK4 is the Runge-Kutta routine 

As written the program does not use an adaptive step-size 

The variables are defined as follows: 

N - the number of ODE's in the system 

NMAX - the number of points saved for plotting 

Currently the program saves every tenth point in a 
file called RKDAT.DAT 
H - the desired stepsize 
X,Y,Z - the independent variables 

Currently the program solves for only three 
variables. If a system has more than three 
equations the additional variables must be 
declared . 

F - the ODE' s 

the equations must be specified in the subroutine 
' DERIVS' 

Tim - time stored as a vector 



All other variables are working storage spaces 



The current program solves the duffing's equation 
with parameters: 

SIG, R, B 

Currently, the output matches the HHC data from Hughes 



The Main Program 

PARAMETER ( N*=2 , NMAX = 9990 ) 

IMPLICIT DOUBLE PRECISION( A-H , O-Z ) 

CHARACTER*? DUMMY 
DATA H, A, B/. 1,7. 5, 0.5/ 

DIMENSION XO(N) ,XEND(N) ,F(N) ,XKl(N) ,XK2(N) ,XK3(N),XK4(N),TO(10) 
COMMON /PARS/H,A,B 



TIME-0.0 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

* set the initial conditions 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

DATA XO/0. 1,4.1/ 
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DUM-U . U 

DUMMY*' AaaaaA' 

0PEN(8,FILE*'dl' , STATUS*' old ' ) 

A 

* the lines of text are added to make the data compatible 

* with the program 'CHAOS' 

A 

WRITE ( 8 , ' ( A ) ' ) ' xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ' 

WRITE ( 8 , ' ( A ) ' ) ' xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ' 

WRITE ( 8 , ' ( A ) ' ) ' xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ' 

DO 11 I-1,NMAX 
X-XO( 1 ) 

Y-XO( 2 ) 

TIME=TIME+H 

A A A A A A A A A A A Af h i s Statement makes the data compatible with 'CHAOS'*********** 

WRITE( 8 , 100 ) DUMMY,TIME,X, Y,X,DUM,DUM 

CALL RK4 (XO,XEND, F,N, TIME,XK1 ,XK2 ,XK3 ,XK4 ) 

DO 14 L-1,N 

XO(L)=XEND(L) 

14 CONTINUE 

12 CONTINUE 

11 CONTINUE 

100 FORMAT (A7 , F8 . 3 , El 1 . 4 , 4 El 2 •. 4 ) 

STOP 

END 

SUBROUTINE RK4 ( XO , XEND , F , N , TIME , XKl , XK2 , XK3 , XK4 ) 

* This subroutine performs numerical integration using 4th order 

* Runge-Kutta methods 



IMPLICIT DOUBLE PREC I S ION ( A-H , 0-Z ) 

COMMON /PARS/H,A,B 

DIMENSION F(N) ,XO(N) ,XEND(N) ,XKl(N) ,XK2(N) ,XK3(N) ,XK4(N) 

DERivs returns values for the ODE'S 
Each loop is a step for N equations 



CALL DERIVS(XO,TIME,F,N) 

DO 51 K*1,N 

XKl(K)=H*F(K) 

XO( K)*XO(K )+XKl( K)/2 . ODO 

51 CONTINUE 

CALL DERIVS(XO,TIME, F,N) 
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XK2(K)=H*F(K) 

X0( K)-XO( K)+XK2( K)/2 . ODO 

52 CONTINUE 

CALL DERIVS( XO,TIME, F,N) 

DO 53 K=1,N 

XK3(K)*H*F(K) 

XO( K)=»XO( K)+XK3( K) 

53 CONTINUE 

CALL DERIVS(XO,TIME,F,N) 

DO 54 K=1,N 

XK4(K)=H*F(K) 

XEND( K)=XEND( K)+(XKl(K)+2. 0D0*XK2( K) +2 . ODO * XK 3 ( K ) +XK 4 ( K ) )/6 . ODO 

54 continue 

RETURN 

END 



SUBROUTINE DER I VS ( XO , T , F , N ) 

This subroutine obtains values for the ODE'S 

This subroutine must be changed for a new system of equations 
This routine solves Dufflng's equation 



IMPLICIT DOUBLE PRECI SION ( A-H , 0-Z ) 
COMMON /PARS/H,A,B 
DIMENSION F(N) ,XO(N) 

F(1)-X0(2) 

F( 2 )-A*COS( T)-XO( 1 ) **3-B*XO( 2 ) 

RETURN 

END 



* use this routine to solve the Lorenz equations 



kAAAAAAAAAAAAA 



bAAAAAAAAAAAAAAAA 



IMPLICIT DOUBLE PRECI S ION ( A-H , 0-Z ) 
COMMON /PARS/SIG,B,R,H 
DIMENSION F(N),XO(N) 



F( 1 )«SIG* ( XO( 2 )-XO( 1 ) ) 

F( 2 )=R*XO( 1 )-XO( 2 )-XO( 1 ) *XO( 3 ) 
F( 3 )-XO( 1 ) *XO{ 2 )-B*XO( 3 ) 



A 

A 

A 



RETURN 

END 



tAAAAAAAAA^ 



tAAAAAAAAAA^ 



^AAAAAAAAA^ 
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